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EUGEN J. lONASCU AND RODRIGO A. OBANDO 

Abstract. In this paper we describe a procedure of calculating the number cubes that have 
coordinates in the set {0, 1, n}. We adapt the code that appeared in [11] developed to calculate 
the number of regular tetrahedra with coordinates in the set {0, 1, ...,n}.. The idea is based on 
the theoretical results obtained in [15]. We extend then the sequence A098928 in the Online 
Encyclopedia of Integer Sequences [16) to the first one hundred terms. 



1. INTRODUCTION 

In this paper we continue and, at the same time revise, some of the work begun in the sequence 
of papers P], [Sj-IIS] about equilateral triangles, regular tetrahedra, and other regular polyhedra, 
all having integer coordinates. Very often we will refer to this property by saying that the various 
objects are in 1? . Strictly speaking these geometric objects are defined as being more than the set 
of their vertices that determines them, but for us here, these are just the vertices. So, for instance, 
an equilateral triangle is going to be a set of three points in 1? for which the Euclidean distances 
between every two of these points are the same. The main purpose of the paper is to take a close 
look at the cubes in 1? . One can imagine easily such cubes by taking the faces parallel to the 
planes of coordinates. However, it is less obvious that there exist many more other cubes sitting 
in space as in Figure [1] (a) . As a curiosity, our counting shows that there are precisely 242,483,634 
cubes with vertices in {0, 1, 100}'^ and one non-trivial example of these cubes is given by the 
points 

{[0,56,59], [21,68,3], [24,0,56], [45,12,0], [52,77,83], [73,89,27], [76,21,80], [97,33,24]} := C. 
In |13] we proved the following theorem. 
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This theorem imphes that there is a one-to-two correspondence between the cubes and the 
regular tetrahedrons in . In [nj we have developed a Maple code to compute the number regular 
tetrahedrons in {0, 1,2, ...,n}'^. We will basically use the same idea and introduce some updates 
based on important theoretical observations. The problem of finding the number of cubes in space 
with coordinates in {0, 1, 2, n} has been studied also in [14J. We list here a few more terms in 
the sequence A098928. 
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It is clear that A098928 < ^103158. For n > 4 we have actually a strict inequality, A098928 < 
j4l03158, and this is due to the fact that some of the tetrahedrons inside of the grid {0, 1, 
extend beyond the grid's boundaries to the unique cube containing it. In Figure 2 we have included 
the two graphs of the sequences A098928 and A103158 up to n = 100. 
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Figure 2: Tetrahedra versus cubes 

2. Theoretical background 

Let us review some of the facts that we are using. Regular tetrahedra are going to be obtained 
from equilateral triangles. Equilateral triangle in are obtained in the following way. Given an 
odd integer d there is a very precise number of solutions for the Diophantine equation (see [11]). 



(1) + + = 3(f , with < a < b < c and gcd{a, b,c) = 1 

which is given by 



(2) 

where 



7re(d) 



A(d) +24r2(rf) 
48 



(3) r^id) = { 



' if d is divisible by a prime factor of the form 8s+5 or 8s+7, s > 0, 
1 if d is 3 

where k is the number of distinct prime factors of d 
of d of the form 8s+l, or 8s+3 (s > 0), 



2" < 
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(4) 



A{d) :=8d Yl 1 

p\d,p prime \ 



p 



and (^) is the Legendre symbol. In particular, these precise countings show that the equation 
([!]) has solutions for every odd d > 1. 

We remind the reader that, if p is an odd prime then 



(5) 



P 



ifp = 3 

1 if p = 1 or 7 (mod 12) 
5 or 11 (mod 12) 



- 1 if p 

As an example, for d = 2011 = 251(8) + 3 = 167(12) + 7, which is a prime, will give A{d) = 16080 
and ^2{d) = 48, and so 7re(2011) = ^^^^^^'^^ = 335. There is only one solution, in this case, for 
which the three values of a, b, and c are not all distinct: a = 139 and b = c = 2461. 

For each solution of ([1]) one can find a two integer parameter family of equilateral triangle in , 
with vertices in a plane through any point of integer coordinates (so we can take simply the origin 
O) and normal ^^0^- Such an equilateral triangle, say AOPQ, can be given terms of two vectors 
C and if described by the next formulae. 



(6) 



oP 



m C 



nlf, 0^ = nt -{n- m)lf, with 't = {(1,(1,(2), if = ivi, 



m,m) 



(7) 



rac + dbs 



das — bcr 



Us 



V3 



db{s — 3r) + ac(r + s) 

da{s — 3r) — 6c(r + s) 



2q 



r + s 



where q = a? + b'^ and (r, s) is a suitable solution of 2q = + 3r^ that makes all the numbers in 



d?]) integers. The sides-lengths of AOPQ are equal to dy^2{m'^ — mn + n^). 

One way to give a more precise construction of a good choice of (r, s) is this to compute the 
greatest common divisor, s + i^/Sr, oi A — i\/^B and 2q in the ring Z[i\/3], where A = ac and 
B = bd. 
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Let us observe that + 3B^ = (ac)^ + 3{bdy = (a^ + b'^){c'^ + 6^) which shows that 2q divides 
+ 3i?2 = (^A + i\/^B){A — i\^B). Since 2q = 4{4k + 1) for some integer k, we are thinking of 
4 as (1 + — i\/3), so the prime factors of 2q here are given by 1 + i\/3, 1 — i\/3 and ah the 

others that appear which are either primes of the form 6fc — 1 or of the form 6k + 1. The factors 
of the form 6k — 1 must appear to even power and those of the form 6k + 1 can be decomposed 
into prime factors u + i\'^v and u — iy/Sv (by Euler's Theorem). Each of these factors can be 
found either in A + i\/^B or A — i^/3B. The product of these factors give s + i-v/3. By construction 
A—i\/3B = {s+i\/3r){u+i\/3v) and 2q = {s+i^/3r){s—i\/3r) because of properties of conjugation 
of complex numbers. This imphes that {A — i\/3B){s — i\/3r)) = 2q{{u + iVSv). 

Hence we get the relations 

2q = s'^ + 3r^, As — 3Br = 2qu and Ar + Bs = —2qv. 
These two relations show that Ci and r/i in ([7]) are integers. Lagrange's identity shows that 



g2 Q,2 _|_ ^2 ' 

which in turn gives 

.2 ^ .2 ^ .2 _ r\a^ + b^) + r^c^ + d^s^ _ d^s^ + 3r^) _ . 
C1+C2+C3 - ^2^2 - - -^<^ ■ 

This implies in particular that C2 must be an integer and that | C I = d\/2- 

It is clear that r and s must be either both odd or both even. This implies that 7^3 is an integer. 
Using again Lagrange's identity we get 



2 I 2 



(a2 + b^)[c\r + s)2 + d\s - 3r)2] c\r + s)^ + d^s - 3rf 



V 4(a2 + 62) 

which implies 



2 2 _ (r + g)2(Q2 + 62 + c2) + dHs - 3r)2 _ d^[3ir + s)^ + {s - 3r)2] _ 

^1+^2+^2- 4(a2 + 62) - Tq -^"^ • 

As before, this proves that 7/2 is an integer and |^| = d^/2. 
In order to find the dot product of C j ^ we observe that 

(a2 + 62)[c2(r2 + r-s) + d2(s2 - 3rs)l 

Cm + C2m = 

which implies 
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_^ _ c^(r^ + rs) + - 3rs) + rs _ (P{3r'^ + 3rs + - 3rs) _ 2 

Hence the angle between the vectors C , V is arccosf S'^ ) = 60°. 
Using these relations, we can easily calculate 

|Oi^P = m^l C p — 2 C • ifmn + n'^\'rf\'^ = 2(f'{rn?' — mn + n^), and 



|0(^P = n^l ^ 2 -if n{n—m)+{n—mf'\l^ ^ = 2d^[n^— m)+(n— m)^] = 2d^ (jm? —mn+n^) . 
The dot product of Oi^ and is then equal to 

Oi^ • = mn\ C P — [m(n — m) + v?] C ■ + — m)\~^\'^ = 

d?{2mn — mn + — + 2n^ — 2mn) = (P{m'^ — mn + n^). 
These relations show that the triangle AOPQ is indeed equilateral and its side lengths are equal 



to — mn + n^. One can easily check that 

aCi + + cCs = aCi + 6??2 + cr/s = 0, 
which implies that AOPQ is indeed contained in the plane of normal it = -^^J^- 

One good question here is, whether or not there are other equilateral triangles with integer 
coordinates contained in the same plane. 

We have proved this in [2j. However, for completion we will include a new relatively simpler 
argument and with a more geometric flavor. If there exists one such triangle, say AOAB, we may 
assume one of its vertices is at the origin (by implementing a translation with integer coordinates). 
Because the vectors C and if form a basis for the space of vectors perpendicular to it, the equation 
OA = m C, —nlf can be solved uniquely for real numbers m and n. Let us consider then the vectors 

OB' = n Q — (n — m)lf, OB^' = [m — n) C, — mlf . 

With the same computations as before this will give two equilateral triangles in the same plane, 
namely AOAB' and AOAB" . Because there are only two equilateral triangles sharing the side 
OA in the given plane, we must have either B' = B 01 B" = B. Without loss of generality, let us 
assume that B' = B. From the formulae in ([7]) we get that 

r+s r+s r+s r—s 

mr n = n G Z, ana nr (n — m) = m H n = u G Z. 

2 ' 2 ^ ^ 2 2 

If we look at these two equations a system of equations in m and n, we get by Cramer's formula, a 
unique solution which are going to be (rational numbers) fractions with integer numerator as the 
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denominators equal to 

r — s r + sr + s + 3r^ q + 5^ 

r 1 = = — = > 1. 

2 22 4 2 2 - 



The same calculations as before show that the side length of IXOAB is given by £ = d\J2[rrfl 
We can do this whole construction for h and c, or for a and c, instead of o and 6. We get that 



(8) = 2S{rr? - mn + r?) = 2cf{mj - mini + nf) = 2(f{ml - m2n2 + nj) 

with mi, ni, m2, n2 rational numbers with ^-^^ and ^-5^ at the denominators, respectively. 
But the gcd( ° ^ , ; " ^'^ ) cannot be greater than one since gcd{a,b,c) = 1. Hence, in ([8]) 
we cannot have the number — mn + = mf — mini -\- nf = m\ — m2n2 + be a fraction 
in the reduced form with a denominator greater then one since any prime dividing it will divide 

2 I-.2 I-.2 2 2 2 

So, we proved that the triangle AOAB (or any other triangle with 
integer coordinates in the plane containing the origin of normal it) has sides at least dy/2. 

Now, if m and n are not integers, then A and B fall strictly inside of the tessellation with 
equilateral triangles generated by the two vectors ( and if (see Figure 3) . Because the tessellation 
is invariant to 60° the position of B inside one of the equilateral triangles is perfectly the same as 
the position of A inside of its equilateral triangle. This creates two vectors of the length and one is 
the rotation of the other by 60°. Using translations with integer coordinates the two vectors show 
the existence of a triangle with one vertex the origin, AOCD, which is equilateral having integer 
coordinates and side lengths strictly less than d\/2. This contradiction shows that A and B must 
be vertices of the tessellation generated by C aiid and so the parametrization (j7]) is unique. 




Figure 3: Two distinct tessellations 



We have given then another proof of Theorem 1 in [lOj . 
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We are including here an example that is illustrating this parametrization. For d = 2011 we have 
seen a solution of ([T]): a = 139 and b = c = 2461. If we do the parametrization with g = + 6^ = 
(2)(3037921) (3037921 is prime). Since A = ac = (23) (107) (139) and B = (23)(107)(2011), we get 



A-iB = (23) (107) (1 - \/3i)(1543 - 468\/3i) 



and 



2q = {l- \/3f)(1543 - 468\/3f)(l + \/3f)(1543 + mVSi) 

which gives s + rVSi = (1 - \/3i)(1543 - mVSi) = 139 - 2011 VSi 
This gives 



Ci = o 

C2 = 2011 and <( 
Cs = -2011 



r]i = -2461 
C2 = 1075 
C3 = -936. 



One can check that, in fact, in the case 6 = c we may always take 



Ci = o 

(9) {^(2 = d and 

= -d 

Let us summarize all these facts that we have shown so far. 



m = -0 

V3 = ^- 



Theorem 2.1. For every solution of the Diophantine equation a^ + b"^ + = 3(i^ there 
exists essentially only one parametrization for the equilateral triangles with vertices in the plane 
of equation ax + by + cz = 0. Every equilateral triangle in 1? is given by such a parametrization. 
Moreover, up to a translation, the vertices of such a triangle are given by (0) and where r and 
s can be computed by finding the greatest common divisor, s + i\/3r, of A — i\/^B and 2q in the 
ring Z[i\/3], where A = ac and B = bd. In the case b = c, the formulae ^ simplify to 

In |12] . we have shown that the only equilateral triangles, in Z^, which can be completed to a 
regular tetrahedron in Z^, are the ones (given as in ([6|) and d?])) for which 'm? — mn + = k"^ for 
some A: € Z. More precisely, if k is divisible by 3 then one can accomplish this on either side of the 
plane containing the triangle and if k is not divisible by 3 then this can be done on only one side. 
By the way, this is saying in particular that, there are a lot more equilateral triangles than regular 
tetrahedrons in 1? . 

The coordinates for the fourth vertex, assuming the equilateral triangle's vertices are as in ([6]) 
and (I?]), are given by 



COUNTING ALL CUBES IN {0,1,..., n}^ 



9 



(2Ci-??i)"i {2C2-V2)m (2C3 - \ 

-{Ci + m)n -iC2 + m)n -{C3 + m)n 

±2ak ±2bk ±2ck 

3 ' 3 ' 3 

/ 

As we already mentioned in Theorem ll.il as long as the coordinates in (|10p are integers then the 
tetrahedron can be completed to a cube in Z^. We are using this formula mostly for k = 1 (let us 
choose m = 1 and n = 0) although there is a need for the general case for big values of d because, 
as pointed out in [IT], there are irreducible regular tetrahedra which cannot be constructed from 
a face as above, by simply taking k = 1. However for small i {£ < 5187 = 3(7)(13)(19)) one can 
find a face of a given regular tetrahedron of sides equal to £\/2 which has the corresponding k as 
in (fTOl) equal to 1. 




/ 

/ 



Figure 4: Eight tetrahedra and essentially one cube 

If one takes all the possible values for m and n such that m? — mn + = 1 there are six regular 
tetrahedra generated this way, from a plane (colored blue in Figure 4), three on one side and the 
other three on the other side, but if one looks at the Figure 4, he/she might observe that in fact 
there are eight regular tetrahedra all generating essentially the same cube (up to translations of 
integer coordinates). As result, our code is going to reflect this property and we will only use one 
of the choices for the values for m and n. Summarizing, there are in general four planes containing 
the center of a given cube in Z^, corresponding to normals given by the directions of the four big 
diagonals in the cube which may generate the the cube as before, some may have a value of > 1. 
For this reason, one needs to check for repetitions when writing the code. For this purpose, our 



(10) 
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approach is to generate an exhaustive hst, £, of cubes in Ng (No = N U {0}) which are irreducible 
(cannot be integer dilated to a smaller cube in Z^). One other property of each cube in L is that it 
cannot be translated in the negative direction along any of the axes of coordinates and remain in 
Nq- However, the cubes in £ are not uniquely defined this way, because of the possible symmetries 
involved here. These symmetries are, in general, 48 in number and form a group which can be 
identified with the symmetry group of a regular octahedron (see p^). 

3. The minimal List and other considerations 




Figure 5: First cubes in C 

A dozen cubes (listed in non-decreasing order of their side-lengths) in the list C that we found 
using our code, are included in the table below. The first column represents the side-lengths, the 
second column gives the dimension of the smallest cube Cm '■= [0, m]^ containing the one in column 
three. 



Table 1 



n 


m 


A cube 


k- values 


invariants 


1 


1 


[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0], [1, 0, 1], [1, 1, OJ, [1, 1, 1] 


1 


[1, 1, 0, 0] 


3 


5 


[0, 3, 2], [1, 1, 4], [2, 2, 0], [2, 5, 3], [3, 0, 2], [3, 3, 5], [4, 4, 1), [5, 2, 3] 


1,3 


[4, 4, 0, 0] 


5 


7 


[0, 0, 4], [0, 5, 4], [3, 0, 0], [3, 5, 0], [4, 0, 7], [4, 5, 7], [7, 0, 3J, [7, 5, 3] 


1 


[12, 18, 4, 0] 


7 


11 


[0, 6, 8], [2, 9, 2], [3, 0, 6], 


5, 3, 




6, 8, 11], [8, 11, 5], [9, 2, 9], [11, 5, 3] 


1,7 


[8, 8, 0, 0] 


9 


15 


[0, 5, 5], [4, 4, 13], [4, 13, 4], 


7, 1, 1 




8, 12, 12], [11, 0, 9], [11, 9, 0], [15, 8, 8] 


1,3 




24, 108, 48, le 




11 


19 


[0, 11, 13], [2, 2, 7], [6, 17, 6], [8, 8, 0], [9, 9, 19], [11, 0, 13], [15, 15, 12], [17, 6, 6] 


1 




24, 108, 48, It 




13 


19 


[0, 12, 15], [3, 16, 3], [4, 0, 12], ]7, 4, 0], [12, 15, 19], [15, 19, 7], 116, 3, 16], [19, 7, 4] 


1, 13 


[8, 8, 0, 0] 


13 


17 


[0, 0, 12], ]0, 13, 12], ]5, 0, 0], [5, 13, 0], ]12, 0, 17], [12, 13, 17], [17, 0, 5], [17, 13, 5] 


1 


[12, 30, 8, 0] 


15 


25 


[0, 5, 10], [2, 19, 15], [10, 0, 20], [11, 7, 0], [12, 14, 25], [13, 21, 5], [21, 2, 10], [23, 16, 15] 


1, 3 


[48, 360, 176, 64] 


17 


29 


[0, 20, 9], [1, 8, 21], [12, 12, 0], [12, 29, 17], [13, 0, 12], [13, 17, 29], [24, 21, 8], [25, 9, 20] 


1 




24, 60, 16, 




17 


23 


[0, 0, 15], [0, 17, 15], [8, 0, 




8, 17, 0], [15, 0, 23], [15, 17, 23], [23, 0, 8], [23, 17, 8] 


1 




12, 42, 12, 




19 


31 


[0, 16, 10], [6, 6, 25], [10, 31, 16 




15, 10, 0], [16, 21, 31], [21, 0, 15], [25, 25, 6], [31, 15, 21] 


1, 19 


[8, 8, 0, 0] 



In the column four we list the values of k which can be used in the construction described in Section 
2 to generate the cube in column three. The list of invariants are as follows. First, is the number 
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of cubes in the orbit obtained by applying the group of 48 transformations determined by the 
orthogonal matrices with coefficients and 1. Let us denote this number by uq. Notice that this is 
a divisor of 48 as expected (Lagrange's Theorem). We expect that in general such a cube will have 
no special symmetry and so, more often we will get = 48. The second number in the invariants 
list is the number of cubes in the generalized orbit, obtained by the previous orbit together with 
all its integer translations along the axes of coordinates that remain in Cm, a number that we are 
going to denote by a. The third number in the list, /3, is the cardinality of the intersection between 
this former orbit and its translation along (0,0,1). Finally, the last number, 7, is defined by the 
cardinality of the generalized orbit with its translation along (0, —1, 1). It turns out that these last 
three numbers are enough to determine the number of cubes that one can fit by translating the 
given cube in all possible ways within a bigger cube of size k > m. This fact has been essentially 
proved in Theorem 2.2 in [lOj. The formula that gives this number is 



(11) {k-m+ if a -3{k-m){k-m+ 1)^/3 + 3(A; - m + 1)(A; - mf'j. 

One of the observations that we will make, about Table 1, is that this set of invariants is not 
complete, since we see that the same numbers appear for various irreducible cubes. The most 
surprising are those cubes in rows six and seven. A good problem here is to determine the exact 
number of such cubes, which go into a certain side-length n, in terms of n. We see that the first n 
for which we have two such cubes is n = 13. Let us also observe that in column four we see a 1 in 
there for each cube. As we mentioned earlier, this is not always the case. 

Also, each cube in Table 1 with side-lengths n, gives rise to an orthogonal matrix with rational 
coefficients having denominators in (obtained by taking the normalized vectors along the sides 
of the cube that form an orthogonal basis). In [I3j we computed a few of them which are included 
here next: 



n ■■= - 



1 -2 

2 -1 
-2 -2 



3 
-4 
-5 



,T7 := - 



-2 -3 
3 -6 
-6 -2 



,Tg:= - 
' 9 



7 


-4 


2 


4 


1 


8 


4 


8 


1 



- IT 



2 -9 -6 
9-2 6 
'6 -6 7 



,ri3 :-- 



-4 -12 -3 
12 -3 -4 
3 -4 12 



,Ti3 :-- 



-13 
12 5 
-5 12 



The next matrix can be obtained by multiplying T3 with T5. We notice a multiplicative structure 
on this set of matrices. For the next two prime sizes we have essentially two orthogonal matrices. 
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15 8 
8 -15 
-17 



15 -6 -10 
10 15 6 . 
6 -10 15 

From T5, T13, and T19 it is clear that there is a natural imbedding of the primitive Pythagorean 
Triples into this sequence of orthogonal matrices (well known in the literature). One interesting 
question here is the following: what is the algebraic relevance for the geometric invariants ao, 
a, (3 and 7, for an orthogonal matrix as above? One more observation here: each cube in the 
list generates by translations and rotations cubes in and two different cubes cannot generate 
the same cube because by doing those transformations the four planes given by the diagonals are 
preserved and those are different for two different cubes in the list. So, in order to count all the 
cubes in C^, we first compute the list of irreducible cubes in £, up to the side length fc, and then 
for each one we use the formula (jlip to find out how many are generated by each in C^. This is not 
enough though because there are cubes in which are not irreducible. So, in the end we multiply 
each cube in the list C by an integer factor in such a way the resulting cube can stil fit in C^. Then, 
we recalculate the invariants on this cube and use the same formula (jlip to find the contribution of 
the reducible cubes. In the end, we add up all these numbers and that gives, NC{k), the number 
of cubes in . 

The first 100 values of iVC are: 1, 9, 36, 100, 229, 473, 910, 1648, 2795, 4469, 6818, 10032, 14315, 
19907, 27190, 36502, 48233, 62803, 80736, 102550, 128847, 160271, 197516, 241314, 292737, 352591, 
421764, 501204, 592257, 696281, 814450, 948112, 1098607, 1267367, 1456292, 1666998, 1901633, 
2162179, 2450440, 2768346, 3117935, 3501389, 3923178, 4384792, 4889323, 5439155, 6037660, 
6687358, 7391669, 8154671, 8979750, 9870158, 10830095, 11862711, 12972046, 14161848, 15436931, 
16801993, 18263634, 19825948, 21493019, 23269647, 25160816, 27171482, 29308957, 31577319, 
33986616, 36540004, 39244371, 42106267, 45131996, 48327502, 51700279, 55258019, 59011634, 
62965766, 67132037, 71515527, 76127374, 80973598, 86062187, 91401297, 96999986, 102866282, 
109014085, 115457359, 122206348, 129266410, 136648555, 144364071, 152426724, 160843660, 169626467, 
178787563, 188347314, 198309846, 208694461, 219509943, 230767760, and 242483634. 

4. The code 

The code is written in Maple and it is attached in pdf format at the end of the paper after 
bibliography. 



Ti 



7 ■- 



1 

17 



12 
12 
1 



9 
-12 



12 



1 

17 



Tig :- 



19 



6 

17 

-6 



-18 1 
6 6 
-1 18 



,^19 :— 



19 
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> restart : with (numtheory) : with (plots) 



1 (Step 1) Procedures to find the possible values 
of k, m and n, and normal plane vector (a,b,c) 

Find all values of k less than n such that k is of the form a product 
of primes of the form 3k+l 

> kvalues : =proc (n) 

> local i, j ,k,L,a,p,q,r,m,inin; 

> L:={};mm:=floor((n+l)/2) ; 

> for i from 2 to mm do 

> a:=ifactors(2*i-l) ; 

> k:=nops(a[2] ) ;r:=0; 

> for j from 1 to k do 

> m:=a[2] [j] [1] ; p:=m mod 3; 

> if m=3 then r:=l fi; 

> if r=0 and p=2 then r:=l fi; 

> od; 

> if r=0 then L:=L union {2*i-l};fi; 

> od; 

> L:=L union {1}; 

> L:=convert(L,list) ; 

> end: 

Find solutions (m,n) of the equation k'2=m'2-mn+n'2, gcd(ni,n)=l,ni>0,n>0,2m<n. 

> listofmn:=proc(k) 

> local a,b,i,nx,x,m,n,L,bl,b2; 

> x:=[isolve(k"2=m"2-m*n+n"2)] ; 

> nx : =nops (x) ; 

> L:={}; 

> for i from 1 to nx do 

> if lhs(x[i] [1] )=m then a:=rhs(x[i] [1] ) ; b : =rhs (x [i] [2] ) ; 

> else b:=rhs(x[i] [1] ) ; a:=rhs(x[i] [2] ) ; fi; 

> if gcd(a,b)=l and a>=0 and b>0 and 2*a<b then L:=L union 

> {[a,b]};fi; 

> od;L; 

> end: 



Find Only the solutions that satisfy gcd(a,b,c)=l and 0<a<=b<=c are in the 
input. 



> abcsol:=proc(d) 

> local i, j ,k,in,u,x,y,sol,cd;sol:={}; 

> for i from 1 to d do 

> u:=[isolve(3*d"2-i~2=x~2+y~2)] ;k:=nops(u) ; 

> for j from 1 to k do 

> if rlis(u[j] [l])>=i and rhs(u[j] [2])>=i then 

> cd:=gcd(gcd(i,rhs(u[j] [1] )) ,rhs(u[j] [2])); 

> if cd=l then sol:=sol union 

> {sort([i,rhs(u[j] [1] ) ,rhs(u[j] [2])])};fi; 

> fi; 

> od; 

> od; 

> convert (sol, list) ; 

> end: 



2 (Step2) The new method of finding r and s. 
It is followed by the implementation of the 
parametrization of equilateral triangles con- 
structing only one regular tetrahedron and 
completing it to a cube (the origin is one of 
its vertices) 

Procedure to get the unique writing as x^2+3y^2 for a prime p of 
the form 6k+l. 

> uniquedecomposition:=proc(p) 

> local out , s , outl , out2 ; 

> if p=2 then out: =[1,1]; fi; 

> if p>2 then s : =isolve (p=x~2+3*y~2) ; 

> outl:=abs(rhs(s[l] [1] ) ) ; out2 :=abs(rhs(s [1] [2])) ; 

> if outl"2+3*out2"2=p then out:=[outl,out2] ; else 

> out:=[out2,outl] ;fi; 

> fi; 

> out ; end : 



The neod procedure is calculating the factorization in Z[sqrt(3)i] of 
a number of the form u+v*sqrt(3)i 



> f actoroverEisensteinintegers : =proc (u, v) 

> local i,M,M,k,a,x,f l,f2,g,y,yl,y2,L,NN,MM,uu,vv; 

> a:=sqrt(3)*I; 

> NN : =gcd (u , v) ; uu : =u/NN ; vv : =v/NN ; 

> N:=uu"2+3*vv~2; 

> x:=uu+vv*a; 

> M: =if actors (N) ;k:=nops(M[2]) ; 

> for i from 1 to k do 

> f 1:=M[2] [i] [1] ; f2:=M[2] [i] [2] ; 

> if fl>2 then 

> g:=uniquedecomposition(f 1) ; 

> else g:=[l,l] ; f2:=l; 

> fi; 

> y:=expand (rationalize (x/(g[l]+a*g [2]))) ; 

> yl:=Re(y); 

> y2 :=type(yl, integer) ; 

> if y2=true then L[i] :=[g[l]+a*g[2] ,f2] ; else 

> L[i]: = [g[l]-a*g[2],f2]; 

> fi; 

> od; 

> [NN , seq(L [i] , i=l . . k) , expand (NM*product (L [ii] [1] "L [ii] [2] , ii=l . . k) ) ] ; 

> end: 

Finding the gcd between A+isqrt(3)B and 2(q) 

> f indgcd:=proc(A,B,q) 

> local 

> i, j ,f .common, m,qq,fac,a,L, Iinportant,nfac,f l,ffl,f 2, s,LL,n,P,rs,LLL; 

> a:=sqrt(3)*I;L:={}; 

> f :=f actoroverEisensteinintegers (A,B) ;m:=nops(f )-l ; 

> common : =gcd (f [1] ,2*q) ; 

> qq:=2*q/common''2; 

> P:=common; 

> f ac : =if actors (qq) ;nf ac : =nops(f ac [2] ) ; 

> for i from 1 to nfac do 

> fl:=fac[2] [i] [1] ;f2:=fac[2] [i] [2] ; 

> if fl=2 then f2:=l;fi; 

> s :=uniquedecomposition(f 1) ; 

> for j from 1 to m do 

> if s[l]+a*s[2]=f [j] [1] then 

> P:=P*(s[l]+a*s[2])-(min(f [j] [2] ,f2)); fi; 

> if s[l]-a*s[2]=f [j] [1] then 

> P:=P*(s[l]-a*s[2])-(min(f [j] [2] ,f2)); fi; 

> od; 

> od; 

> P:=expand(P) ; 

> rs:=[Re(P) ,Im(P)/sqrt(3)] ; 

> [rs,A*rs[l]+3*B*rs[2] mod 2*q, A*rs [2] -B*rs [1] mod 2*q] 

> end: 

The fourth point to get a regular tetrahedron is the origin. We have included 
only one since all the other are going to produce the same cube when transformed 
into the positive quadrant. 



> f indpar:=proc(a,b,c,inm,nn) 

> local i , j , sol , mx , nx , r , s ,my , ny , q , d , u , V , w , X , y , z , mu , nu , mv , nv , ef , ns , mz , nz 

> ,mw,nw,T,l,uu,t,Rl,R2,fc,k,A,B,rs,G,DD,E,F; 

> q:=a"2+b~2;k:=sqrt(inin~2-inin*nn+nn~2) ;d:=sqrt((a~2+b~2+c~2)/3) ;A:=a*c;B: 

> =b*d; 

> rs:=f indgcd(A,-B,q) ; 

> r:=rs[l] [2] ; s:=rs[l] [1] ; 

> uu:=(s"2+3*r"2-2*q)"2; if uu>0 then t:=s;s:=r;r:=t; fi; 

> mx:=-(d*b*(3*r+s)+a*c*(r-s))/(2*q) ;nx:=-(r*a*c+d*b*s)/q; 

> my:=(d*a*(3*r+s)-b*c*(r-s))/(2*q) ;ny:=-(r*b*c-d*a*s)/q; 

> mz:=(r-s)/2;nz:=r; 

> mu : =nx ; mv : =ny ; mw : =nz ; nu : =nx-mx ; nv : =ny-my ; nw : =nz-inz ; 

> u : =mu*m-nu*n ; v : =mv*m-nv*n ; w : =mw*m-nw*n ; 

> x : =mx*m-nx*n ; y : =my*m-ny*n ; z : =mz*m-nz*n ; 

> R2:=[(x+u-2*a*k)/3, (v+y-2*b*k)/3, (z+w-2*c*k)/3] ; 

> Rl:=[(x+u+2*a*k)/3, (v+y+2*b*k)/3, (z+w+2*c*k)/3] ; 

> f c :=subs(m=mm,n=nn,Rl [1] ) ; 

> if fc=floor(fc) then 

> T:=subs(m=mm,n=nn, [[u,v,w] , [x,y,z] ,R1] ) ; else 

> T:=subs(m=mm,n=nn, [[u,v,w] , [x,y,z] ,R2] ) ; 

> fi; 

> G: = [(T[1] [1]+T[2] [1]+T[3] [1])/2,(T[1] [2]+T[2] [2]+T[3] [2])/2,(T[l] [3] + 

> T[2] [3]+T[3] [3])/2] ; 

> DD: = [(T[1] [1]+T[2] [1]-T[3] [1])/2,(T[1] [2]+T[2] [2]-T[3] [2])/2,(T[l] [3] 

> +T[2] [3]-T[3] [3])/2] ; 

> E : = [ (T [1] [1] -T [2] [1] +T [3] [1] ) /2 , (T [1] [2] -T [2] [2] +T [3] [2] ) /2 , (T [1] [3] - 

> T[2] [3]+T[3] [3])/2] ; 

> F: = [(-T[l] [1]+T[2] [1]+T[3] [1] )/2, (-T[l] [2]+T[2] [2]+T[3] [2] )/2, (-T[l] [ 

> 3]+T[2] [3]+T[3] [3])/2]; 

> [T[l] ,T[2] ,T[3] ,G,DD,E,F] ; 

> end: 



3 (Step 3) Procedures for adding vectors, other 
operations on vertices, finding the four k- values 
of a cube (most of the time is less than four) 

> subtrv : =proc (U, V) 

> local W; 

> W[l] :=U[1]-V[1] ;W[2] :=U[2]-V[2] ;W[3] :=U[3]-V[3] ; 

> [W[l] ,W[2] ,W[3]] ; 

> end: 

> addvec:=proc(U,V) 

> local W; 

> W[l] :=U[1]+V[1] ;W[2] :=U[2]+V[2] ;W[3] :=U[3]+V[3] ; 

> [W[l] ,W[2] ,W[3]] ; 

> end: 



> addvectCube:=proc(T,v) 

> local i,Q;Q:={}; 

> for i from 1 to 8 do 

> Q:=Q union {addvec(T[i] ,v)}; 

> od; 

> Q; 

> end: 

> distance :=proc (A, B) 

> local C; 

> C:=subtrv(A,B) ; 

> sqrt(C[l]-2+C[2]-2+C[3]~2); 

> end: 

> f inddiagonall : =proc (C) 

> local i ,m,L,x,y,CC; 

> L:=convert(C,list) ;in:=nops(L) ; 

> for i from 1 to m do 

> x:=type (distance (L[l] ,L[i] )/sqrt(3) .integer) ; 

> if x=true then y:=[L[l] ,L[i]] ;f i; 

> od; 

> CC:=C minus {y[l],y[2]}; 

> y,CC; 

> end: 

> f inddiagonal:=proc(C) 

> local x,xl,x2; 

> x:=f inddiagonall(C) ; 

> xl :=f inddiagonall (x [2] ) ; 

> x2 : =f inddiagonall (xl [2] ) ; 

> [x[l] ,xl[l] ,x2[l] , convert (x2 [2] ,list)] ; 

> end: 

> unitvector:=proc(U) 

> local i, j ,k,l,x; 

> i:=U[l];j:=U[2];k:=U[3]; 

> l:=gcd(gcd(i, j) ,k) ; 

> x:=(i"2+j"2+k"2)/(3*l-2) ; 

> sqrt (x) ; 

> end: 

> f ourkvalues : =proc (C) 

> local Nl, N2,N3,N4,x, length; length : =distance (C [1] ,C [2]) ; 

> X : =f inddiagonal (C) ; 

> Nl:=unitvector(subtrv(x[l] [l],x[l] [2])); 

> N2:=unitvector(subtrv(x[2] [1] ,x[2] [2])); 

> M3:=unitvector(subtrv(x[3] [l],x[3] [2])); 

> M4:=unitvector(subtrv(x[4] [1] ,x[4] [2])); 

> {length/Ml , length/N2 , length/N3 , length/N4} ; 

> end: 

> multbyf actorv:=proc(v,k) 

> local w; 

> w[l] :=v[l]*k;w[2] : =v [2] *k ; w [3] :=v[3]*k; 

> [w[l] ,w[2] ,w[3]] ; 

> end: 



> multbyf actor :=proc(T,k) 

> local i,MT,Q;MT:={}; 

> Q : =convert (T , list) ; 

> for i from 1 to 8 do 

> NT:=NT union {multbyfactorvCQCi] ,k)}; 

> od;NT; 

> end: 

4 (Step 4) Translation of the cube obtained in 
Step 2 into the first quadrant and getting the 
orbit and the generalized orbit, then intersec- 
tions with its translations along (0,0,1) and 
(0,1,-1). Calculation of alpha_0, alpha, beta 
and gamma for a cube. 

Translating the cube to the first octant (x>—0,y>—0,z>=0) 

> tmttopqcube:=proc(T) 

> local i,a,b,c,v,0,C; 

> a:=min(T[l] [1] ,T[2] [1] ,T[3] [1] ,T[4] [1] ,T[5] [1] ,T[6] [1] ,T[7] [1] ,0) ; 

> b:=min(T[l] [2] ,T[2] [2] ,T[3] [2] ,T[4] [2] ,T[5] [2] ,T[6] [2] ,T[7] [2] ,0) ; 

> c:=min(T[l] [3] ,T[2] [3] ,T[3] [3] ,T[4] [3] ,T[5] [3] ,T[6] [3] ,T[7] [3] ,0) ; 

> 0:=[0,0,0]; 

> v:=[a,b,c] ;C:= 

> {subtrv(0,v) ,subtrv(T[l] ,v) ,subtrv(T[2] ,v) ,subtrv(T[3] ,v) , 

> subtrv(T[4] , v) , subtrv(T [5] ,v) , subtrv(T [6] , v) , subtrv(T [7] ,v)}; 

> end: 

> mscofmcube:=proc(Q) 

> local a,b,c,T; 

> T : =convert (Q , list) ; 

> a:=max(T[l] [1] ,T[2] [1] ,T[3] [1] ,T[4] [1] ,T[5] [1] ,T[6] [1] ,T[7] [1] ,T[8] [1] 

> ); 

> b : =max (T [1] [2] , T [2] [2] , T [3] [2] , T [4] [2] , T [5] [2] , T [6] [2] , T [7] [2] , T [8] [2] 

> ); 

> c : =max (T [1] [3] , T [2] [3] , T [3] [3] , T [4] [3] , T [5] [3] , T [6] [3] , T [7] [3] , T [8] [3] 

> ); 

> max(a,b,c); 

> end: 



orbit ICube : =proc (T) 
local 

i,k,Tl,a,b,c,x,T2,T3,T4,T5,T6,T7,T8,T9,T10,Tll,T12,T13,T14,T15,T16,T17 

,T18, 

T19,T20,T21,T22,T23,T24,S,Q,d,al,bl,cl; 
Q:=convert(T,list) ; 
d : =mscof mcube (T) ; 
T1:=T; 

T2:={seq([Q[k] [2] ,Q[k] [3] ,Q[k] [1]] ,k=l. .8)}; 
T3:={seq([Q[k] [1] ,Q[k] [3] ,Q[k] [2]] ,k=l . .8)}; 
T4:={seq([q[k] [1] ,Q[k] [2] ,d-Q[k] [3] ] ,k=l . . 8) } ; 
T5:={seq([Q[k] [2] ,Q[k] [3] ,d-Q[k] [1] ] ,k=l . . 8) } ; 
T6:={seq([Q[k] [1] ,Q[k] [3] ,d-Q[k] [2] ] ,k=l . . 8) } ; 
T7:={seq([Q[k] [1] ,d-Q[k] [2] ,Q[k] [3] ] ,k=l . . 8) } ; 
T8:={seq([Q[k] [2] ,d-q[k] [3] ,Q[k] [1] ] ,k=l . . 8) } ; 
T9:={seq([Q[k] [1] ,d-q[k] [3] ,Q[k] [2] ] ,k=l . . 8) } ; 
T10:={seq([d-Q[k] [1] ,C![k] [2] .QLk] [3]] ,k=l. .8)}; 
Tll:={seq([d-Q[k] [2] ,Q[k] [3] ,Q[k] [1]] ,k=l. .8)1; 
T12:={seq([d-q[k] [1] ,q[k] [3] ,q[k] [2]] ,k=l. .8)}; 
T13:={seq([q[k] [1] ,d-q[k] [2] ,d-q[k] [3]] ,k=l. .8)}; 
T14:={seq([q[k] [2] ,d-q[k] [3] ,d-q[k] [1]] ,k=l. .8)}; 
T15:={seq([q[k] [1] ,d-q[k] [3] ,d-q[k] [2]] ,k=l. .8)}; 
T16:={seq([d-q[k] [1] ,d-q[k] [2] ,q[k] [3]] ,k=l. .8)}; 
T17:={seq([d-Q[k] [2] ,d-q[k] [3] ,q[k] [1] ] ,k=l . . 8) } ; 
T18:={seq([d-Q[k] [1] ,d-Q[k] [3] ,Q[k] [2] ] ,k=l . . 8) } ; 
T19:={seq([d-q[k] [1] ,q[k] [2] ,d-q[k] [3]] ,k=l. .8)}; 
T20:={seq([d-q[k] [2] ,q[k] [3] ,d-q[k] [1]] ,k=l. .8)}; 
T21:={seq([d-q[k] [1] ,q[k] [3] ,d-q[k] [2]] ,k=l. .8)}; 
T22:={seq([d-q[k] [1] ,d-q[k] [2] ,d-q[k] [3] ] ,k=l . . 8) } ; 
T23:={seq([d-q[k] [2] ,d-q[k] [3] ,d-q[k] [1] ] ,k=l . . 8) } ; 
T24:={seq([d-Q[k] [1] ,d-q[k] [3] ,d-q[k] [2]] ,k=l. .8)}; 
S: = 

{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17,T18,T19,T 

20,T21,T22,T23,T24}; 
S; 

end: 

orbitCube : =proc (T) 
local S,q,Tl; 
q:=convert(T,list) ; 

Tl : ={seq( [Q [k] [3] , q [k] [2] , q [k] [1] ] ,k=l . . 8) } ; 
S:=orbitlCube(T) union orbitlCube(Tl) ; 

S; 

end: 



> translCube:=proc(T) 

> local S,Q,i, j ,k,a2,b2,c2,a,b,c,d,alpha; 



> 


Q 


=convert(T,list) ; 




> 


a 


=max(seq(Q [k] [1] ,k=l. 


.8)) 


> 


b 


=max(seq(Q [k] [2] ,k=l. 


.8)) 


> 


c 


=max(seq(Q [k] [3] ,k=l. 


.8)) 


> 


d 


=max(a,b, c) ; 





> a2:=d-a;b2:=d-b;c2:=d-c; 

> S : =orbitCube (T) ; alpha : =nops (S) ; 

> for i from to a2 do 

> for j from to b2 do 

> for k from to c2 do 

> S:=S union orbitCube(addvectCube(T, [i, j ,k] )) ; 

> od; 

> od; 

> od; 

> #alpha,nops(S) , [a2,b2,c2] , (a2*b2+b2*c2+c2*a2)*alpha; 

> S; 

> end: 

> alphabetagamma: =proc (T) 

> local X, S,i, SI, S2,d, cube, j ,exitv, alpha, beta, gamma; 

> d:=mscofmcube(T) ;x:=nops(orbitCube(T)) ; 

> Sl:=translCube(T) ;S2:=Sl;S:=convert(S2,list) ;alpha:=nops(S) ; 

> for i from 1 to alpha do 

> cube:=S[i] ; 

> j:=l; exitv:=0; 

> while j<=8 and exitv=0 do 

> if cube[j] [3] =d then exitv:=l; 

> S2:=S2 minus {cube};fi; 

> j:=j+i; 

> od; 

> od; 

> beta:=nops(S2) ;S:=convert(S2,list) ; 

> for i from 1 to beta do 

> cube : =S [i] ; 

> j : =1 ; exitv : =0 ; 

> while j<=8 and exitv=0 do 

> if cube[j] [2] =d then exitv :=1; 

> S2:=S2 minus {cube};fi; 

> j:=j+i; 

> od; 

> od; 

> gamma: =nops (S2) ; 

> [S2, [x, alpha, beta, gamma] ] ; 

> end: 



5 (Step 5) The function that gives the number of 
cubes generated by a small cube which fits in 
a cube of dimension d, in an arbitrary bigger 
cube of dimension n. 

> f : = (n , d , alpha , beta , gamma) -> (n-d+l) ~3*alpha-3* (n-d) * (-d+l+n) ~ 2*beta+3* 

> gamma* (n-d+ 1 ) * (n-d) ~ 2 : 

6 (Step 6) List of minimal (irreducible) cubes 
together with n, minimal size containg the 
cube, the k- values, alpha, beta, and gamma. 
In the end another list is created in a new pro- 
cedure to get the list of the reducible cubes 
together with the alpha, beta, and gamma. 

We compute the list of minimal cubes taking into account that at 
some point only k=l is not sufficient to generate a cube . 

> ExtendList :=proc(n,N,L: : list, mm, nn, Orb: :array) 

> local i, sol, iisol,nel,ttpcube, orb :: array, 

> LL: :list,tnel,NL: :list,cio,C,pC,k,kvalues,m, exception, abg, MM; 

> NN:=floor((N-l)/2) ; 

> orb:=array(l. .2*NM+1) ; 

> nel : =nops (L) ; 

> LL : =L ; 

> k:=sqrt(mm"2-mm*nn+nn~2) ; 

> m:=n*k; 

> if m<=N then 

> sol:=abcsol(n) ;nsol:=nops(sol) ; 

> tnel:=nel;orb:=Orb; 

> for i from 1 to nsol do 

> C : =f indpar (sol [i] [1] , sol [i] [2] , sol [i] [3] ,mm,nn) ;pC : =tmttopqcube(C) ; 

> kvalues : =fourkvalues (pC) ; exception: =evalb(l in kvalues) ; 

> if k=l or exception=f alse then 

> cio:=evalb(pC in orb[m]); 

> if cio=false then 

> abg:=alphabetagcimma(pC) ; 

> orb [m] : =orb [m] union abg[l]; 

> NL[tnel+l] :=[n,mscofmcube(pC) ,pC, kvalues, abg [2] ,sol[i]] ; 

> tnel : =tnel+l ; 

> fi; 

> fi; 

> od; 

> LL:=[seq(L[i] ,i=l. .nel) ,seq(NL[j] ,j=nel+l. .tnel)] ; 

> fi; 

> LL , orb ; 

> end: 



> ExtendListuptoN:=proc(N) 

> local i, j ,k,l,kv,kvn,NN,L,Orb: :cirray,inn,iiinn,E,n,ii; 

> kv:=kvalues(N) ;kvii:=nops(kv) ; 

> NN:=floor((N-l)/2) ; 

> L: = [] ;Orb:=array(l . .2*NN+1) ; 

> for ii from 1 to 2*NN+1 do 

> Orb[ii]:={}; 

> od; 

> for i from 1 to kvn do 

> k:=kv[i] ;mn:=listofmn(k) ; mnn : =nops (mn) ; 

> for j from 1 to mm do 

> for 1 from 1 to NN+1 do 

> n:=2*l-l; 

> E:=ExtendList(n,N,L,mn[j] [1] ,mn[j] [2] ,Orb) ; 

> L:=E[1]; 

> for ii from 1 to 2*NN+1 do 

> Orb[ii] :=E[2] [ii] ; 

> od; 

> od; 

> od; 

> od; 

> L; 

> end: 

> L:=ExtendListuptoN(50) : 

> ExtendListuptoNmultiples : =proc (N,L) 

> local i, j ,x,lc,m,mm,d,dd,C,CC,LL; 

> m:=nops(L) ;i:=l;LL:={}; 

> while i<=m do 

> d:=L[i][2]; 

> if d<=N then 

> mm : =f loor (N/d) ; C : =L [i] [3] ; j : =2 ; 

> while j<=mm do 

> CC:=multbyf actor (C,j) ;x:=alphabetagamma(CC) ; 

> dd:=d*j ;lc:=nops(LL) ; 

> LL:=LL union { [L [i] [1] *j ,dd,CC,L [i] [4] ,x [2] ] } ; 

> j:=j+i; 

> od; 

> fi; 

> i:=i+l; 

> od; 

> convert (LL , list) ; 

> end: 

> LL:=ExtendListuptoNmultiples(50,L) : 



7 (Step 7) Adding up the contribution of each 
cube. 

> addupnew:=proc(N,L,LL) 

> local i, j ,k,nc,inm,m,d,C,CC,x,dd,nt; 

> nc:=0; 

> ni:=nops(L) ; 

> i:=l; 

> while i<=m do 

> d:=L[i][2]; 

> if d<=N then 

> nc : =nc+f (N , d , L [i] [5] [2] , L [i] [5] [3] , L [i] [5] [4] ) ; 

> fi; i:=i+l; 

> od; 

> m:=nops(LL) ; 

> i:=l; 

> while i<=m do 

> d:=LL[i][2]; 

> if d<=N then 

> nc : =nc+f (N , d , LL [i] [5] [2] , LL [i] [5] [3] , LL [i] [5] [4] ) ; 

> fi; i:=i+l; 

> od; 

> nc; 

> end: 

> NC:=[seq(addupnew(k,L,LL) ,k=l. .50)] ; 

NC := [1, 9, 36, 100, 229, 473, 910, 1648, 2795, 4469, 6818, 10032, 14315, 19907, 
27298, 36886, 49133, 64531, 83784, 107542, 136551, 171599, 213524, 263202, 
321849, 390415, 469932, 561492, 667305, 789317, 929098, 1088500, 1269367, 
1473635, 1703708, 1961706, 2251289, 2575291, 2936272, 3337026, 3780455, 
4269605, 4813854, 5414560, 6076915, 6804587, 7603120, 8476390, 9430481, 
10471175] 



